Fluid Modes of a Spherically Confined Yukawa Plasma 
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The normal modes of a three-dimensional Yukawa plasma in an isotropic, harmonic confinement 
are investigated by solving the linearized cold fiuid equations. The eigenmodes are found analytically 
and expressed in terms of hypergeometric functions. It is found that the mode frequencies solely 
depend on the dimensionless plasma parameter ^ = nR, where R is the plasma radius and k the 
inverse screening length. The eigenfrequencies increase monotonically with ^ and saturate in the 
limit ^ — ^ CX3. Compared with the results in the Coulomb limit [D. H. E. Dubin, Phys. Rev. Lett. 
66, 2076 (1991)], we find a new class of modes characterized by the number n which determines the 
number of radial nodes in the perturbed potential. These modes originate from the degenerate bulk 
modes of the Coulomb system. Analytical formulas for the eigenfrequencies are derived for limiting 



PACS numbers: 52.27.Lw,52.35.Fp,52.27.Gr 



I. INTRODUCTION 

The interaction between charged particles is known to 
be strongly affected by a background plasma. Examples 
include dusty plasmas, where the screening of the dust- 
dust interaction is mainly determined by ions [TH3]; and 
dense two-component plasmas [IHZ], where the ions are 
screened by weakly coupled electrons. These plasmas are 
expected to occur in the interior of giant planets and 
white dwarf stars. While in the former case the degree 
of screening is determined by the ion Debye length, the 
screening length in the latter is the Thomas- Fermi length, 
owing to the degeneracy of the electrons. The Yukawa 
one-component plasma model is often used to describe 
the heavy plasma component while the light component 
determines the screening length. Its static properties and 
collective excitations have been studied in several publi- 
cations, e.g. [STOJ]. 

In many situations the plasma is neither homogeneous 
nor macroscopic. Recently, it was shown [TTl [T^] that 
the density of a three-dimensional dusty plasma, where 
gravity is balanced by a thermophoretic force |13j , is not 
homogeneous. The reason is the screened dust-dust inter- 
action which produces an inhomogeneous density profile 
in a harmonic confinement. This is different from exper- 
iments with confined ions |14j where the interaction is 
Coulombic and the mean density is constant. In astro- 
physical plasmas the confinement is provided by gravity 
and may also influence the plasma properties. 

Previous continuum theories [TTl HI] for Yukawa plas- 
mas were limited to static properties. Here we extend 
these results to a time-dependent theory and investigate 
the normal modes of a Yukawa plasma in a spherical, 
harmonic confinement (T^l [TS] . This model is appropri- 
ate for the experiments of [T31 [T7], for which the nor- 
mal modes of rather small dust crystals have recently 
been measured [18]. On the one hand, a fiuid approach 
is expected to be accurate for long wavelength modes 



in a weakly coupled plasma. On the other hand, the 
agreement of theoretical predictions [19j with experi- 
ments [201 EI] and simulations [22] for confined ions turns 
out to be surprisingly good even in the strongly coupled 
phase. An analogous result for a confined one-component 
plasma with a screened interaction is still missing. Open 
questions are the influence of screening on the normal 
modes and the eigenfrequencies. Compared with the La- 
grangian description of Ref. [23 for the breathing mode, 
the present approach makes no assumption about the 
particular mode form. Besides dusty plasmas and com- 
pact star interiors, we expect our results to be relevant for 
other systems as well, when screening and confinement 
are not negligible. 

This paper is organized as follows. The fluid equa- 
tions are introduced and linearized in Sec. [Ill In Sec. IIIII 
we explicitly consider an isotropic harmonic confinement. 
The density profile is reviewed and used to calculate 
the ground state potential and energy. Further, the lin- 
earized Poisson equation is solved and the eigenfrequency 
spectrum is derived. The normal modes are discussed in 
detail. We conclude with a discussion of the theory and 
an outlook on future work in Sec. IIVI 



II. FLUID DESCRIPTION 

A. Basic equations 

The fiuid equations for a spatially confined one- 
component plasma read 

dn 
iH 
dv 



+ V • (nv) = 0, 



(la) 



dt 



(v- V)v 



-nVC/ — V • P — rrmvw, (lb) 



where f/(r) = V{v) + q(j){r) denotes the sum of the con- 
finement potential T^(r) and the potential (j){r,t) induced 



by the particles. In the first (continuity) equation n{r,t) 
is the particle density and v(r, i) their mean velocity. 
The second equation is the momentum equation, where 
m denotes the particle mass, q their charge and P{r,t) 
the pressure tensor. A damping term with friction coef- 
ficient v is included to account for collisions with neutral 
particles. 

The fluid equations are complemented by Poisson's 
equation for the (induced) potential 0, 



A- 



-Airgn, 



(2) 



where the screening of the interaction between the heavy 
particles by a polarizable background medium (light 
charged components) is explicitly taken into account. 
The range of the interaction is determined by the inverse 
of the screening parameter k. 



where the plasma dielectric function is given by 



eir,io) = l- 






and 



Wp(r) = ^47rg2no(r)/m 



(7) 



(8) 



denotes the local plasma frequency, fii and vi have been 
eliminated in favor of 0i . 

Eq. ([6| is a self-contained equation for (f>i and will be 
solved in the following section for a special case. Having 
found its solution, fii and vi follow from Eqs. (l5|. 



III. SOLUTION FOR HARMONIC 

CONFINEMENT 



Linearization 



Ground state 



A small perturbation of the plasma equilibrium is 
well described by linear response theory. This could 
be caused by an external perturbation (e.g. laser ma- 
nipulation of particles in a dusty plasma) or by ther- 
mal effects. We are interested in strongly coupled plas- 
mas and hence can neglect the pressure term (cold 
fluid limit). Eqs. (TJ2) are then linearized according to 
n(r,t) ~ no{r) +njjr,t), v(r,i) ~ Vi(r,i), (j){r,t) ~ 
4>o{r) + (j)i{r,t). Products of first order terms are as- 
sumed negligible. 

The classical equilibrium density profile noir) and the 
associated potential (j)Q{r) are determined from the zero 
order terms of Eqs. ( Ibj 2]), 



gV0o 



-vy, 



(A - K^) 00 = -47rqno, 



(3a) 
(3b) 



which describe local force equilibrium and are equivalent 
to the energy minimization in |llj (we neglect the finite 
size factor (TV — l)/iV, where N is the particle number). 
First order quantities are determined by 



drii 
~dt 



+ V-(noVi) = 0, 



9vi 
(A-K^ 



1 = —Airqui. 



(4a) 

(4b) 
(4c) 



Looking for normal mode solutions with a time depen- 



dence 



, e.g. (t)i{v,t) — 0i(r)e , we obtain 



iLoni 
m{(jj + ij^)vi 



V • (novi), 
—iqV(j)i. 

Using Eqs. ^ we can rewrite (4c) as 

V- [e(r,cj)V0il =K^^i, 



(5a) 
(5b) 



(6) 



So far our results are valid for arbitrary confinement. 
In order to make further progress let us now explic- 
itly consider an isotropic harmonic confinement V{r) = 
miJ^r^ jl. The ground state density no(r) [cf. Eqs. (|3|] 
is given by [11] 



no(r) 



3 



1 



6 1 



6 



■je(i? 



0, (9) 



where a — ((j^ jmuj^Y^'^ is the Wigner-Seitz radius in the 
Coulomb limit, k = 0. The normalized cluster radius is 
denoted by ^ = kR. For Coulomb interaction the den- 
sity is constant, f = 0, and R{k = 0) = Re = aN^I'^ , 
while for K 7^ no(r) decreases parabolically towards the 
boundary. In this case R(k) = ^/k must be determined 
from [n] 



e6 + 6e' + 15[C^-K^3-4(^ + l)]=0, 



(10) 



where he — uRc is the inverse screening length normal- 
ized by the Coulomb radius. For small fee the asymptotic 
solution of Eq. ([To| is 



£,{kc) ^kc- —K 



15 



9 '^ 25 '^ 



(11) 



while for fcc ^ 1 



i{kc) ^ 15I/5 4/' - 1 + 



1 



152/5 "C7 

25 '^ 15 ^ 



^-6/5 



(12) 



for 



The relative error of these approximations is < 10^^ 
kc < 1.26 and kc > 1.26, respectively. 

The case ^ ^ 1 is encountered if kc — kgN^^^ ^ 1, 
i.e. if either kq and/or N are large. This is why we 
will refer to ^ —^ 00 as the macroscopic/strong screening 
limit. The plasma has a size of many screening lengths. 



The opposite case, ^ ^ 1, will be referred to as the 
Coulomb liniit, where the screening length is much larger 
than the plasma radius. 

It is straightforward to calculate the moments of the 
density, which are given by 



1 

N 



r" no{r) dr 



(13) 



_ IT /Ry e + {n + 6)e + 3(n + 5)(1 + Q 
~ N \a) (n + 5)(n + 3)(l + e) 

The ratio of two moments [Fig. nk] could help determine 
the unknown parameters ^ and R in experiments, where 
the particle positions are directly accessible. The mo- 
ments can easily be calculated since the integral reduces 
to a sum over all particles as n{r) — J2i '^(r — r^). 

The ground state potential (l)o{r) is determined by 
Eq. (3b) for which the Yukawa potential is the associ- 



ated Green's function [TT]. Thus, the solution is given by 
(details can be found in Appendix [A) 



Mr) = q 



-K|r-r'| 

no{T')-. ^dv' 




i[3 + e-(l + 0[i]' 
i?exp(^ — Kr)/r, 



(14) 

r<R, 
r>R. 



Since the confinement is parabolic the potential inside 
the plasma must decrease correspondingly to ensure force 
equilibrium. Outside the cloud the potential behaves like 
that of a point charge placed at the origin. While for 
Coulomb interaction the (effective) charge is Qcff — Nq, 
as expected from Gauss's law, the result for Yukawa in- 
teraction is Qoff — q [R/aY e^/[l + ^]. 

The ground state density and potential can further be 
used to calculate the total energy in mean-field approxi- 
mation, Etot = -Spot + -Bint, where [11] 



-B, 



pot ^ / V{r)na(r)dv, E^nt ^ ^ I (t>o{r)no{r)dr. 



Using Eqs. ( 9|14 | we find after some algebra, 

E, 



^pot ^ {R/af 

q^/a 10 

-Bint ^ jR/a)^ 

q^/a 210 



, e'8 + n 
7 1 + e. 
126 + 147^ + 72^2 ^ 18,^3 



(15) 



2^4 



The total energy then reads 



Et. 



q^ja 



[Rial 
210 



189 -t- 273C -h 159^2 ^ 45^3 _^ 5^4 



(1 + C)^ 



see Fig. [Tja. For small <^ the interaction energy yields 
the dominant contribution to the total energy since the 
potential is only weakly screened. For large ^ the cluster 
has a size of several screening lengths and the potential 
energy dominates. The critical point is at ^ « 1.72. 
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FIG. 1: Screening dependence of a) ratio of n-th to first den- 
sity moment, and b) ground state energy contributions. 



For Coulomb interaction (^ ~ 0) the result is i?tot = 
9/10 iV^/^ (?2/a, which is the first (mean-field) term in the 
energy expansion of the shell model (or the energy of the 
neutralizing background) PHHS] . In the opposite limit, 
f ^ 1, the asymptote is 



E,,. 



R^e 1 e^ 15^/' N^^ 



q^/a \aj 42 (Ka)5 42 42 (Ka)4/5 



, (16) 



since ^ ~ 15"'^/^ k^ [leading order term in Eq. (12)]. In 



the limit iV :^ 1 we thus have Etot oc N^/"^ for k = and 
Etot oc h^^/^N'^/^ for finite screening. 

B. Normal modes 

Since the ground state density profile terminates in a 
finite step, cf. Eq. (|9]), one has to solve ^ separately 
for r < R and r > R. The dielectric function inside and 
outside the plasma reads 



e(Kr, iu) 



where the constant term is given by 




(17) 



ei(w)-l 



= 1- 



^^(0) 



1 



^^(0) 



uj{uj + iv) 

If. e'i 



n^ 



02 



2 1 



(18) 



and fl^ — uj{u! + iv^jiJ^. Accordingly we define a nor- 
malized plasma frequency ^pinr) — u>p{r)/ujQ. Note that 
the plasma frequency ^p{r) only depends on the product 
Acr, see Eqs. (8|9). 

Eq. (|6]) must be supplemented by the boundary condi- 
tions 



efV0i"(r,w) 
Cr ■ e(Kr, a;)V0'"(r,a;) 



r=R 



e.-v^rw 



lout / 



er-v^r'w 



r=R 



r=R 



lim (/)°"*(r)=0. 



(19a) 
(19b) 
(19c) 



Here e* and ir are unit vectors in the tangential and ra- 
dial direction at the surface of the sphere with radius R. 
These are the usual boundary conditions for the tangen- 
tial component of the electric field — V0i and the radial 
component of — eV0i. 

In order to solve Eq. (l7| we use an expansion in spher- 
ical harmonics, i.e. 

Since the spherical harmonics are eigenfunctions of the 
angular part of the Laplacian, 

this leads to the following equation for the radial function 



d_ 
dx 



e{x,uj)x f{x,uj) 



(20) 



[x^ + l{i ^ \)e{x,u)\ f{x,uj) = Q, 



after multiplying by x^ . Here we made a change of vari- 
ables from r to the dimensionless radius x = nr and 
introduced a new notation /(r) -> f{x). In the remain- 
der of this section we will separately solve ( [20| ) inside and 
outside the plasma. 

Consider first the situation outside the plasma where 



the dielectric function is just a constant. Here Eq. (20) 
reduces to 



„2 rll 



f"{x) + 2xf{x) - [x^ + i{t + 1)] /(x) = 0, (21) 

the solutions of which are modified spherical Bessel func- 
tions of the first and second kind, i^(a;) and ke{x), re- 
spectively [26_. They are related to the modified Bessel 
functions by 



ii{x) 



Ie+i/2ix), h{x) 



Kt+i/2{x). 



Only ki{x) is compatible with the boundary condi- 
tion ( 19c I and goes to zero at infinity, so the solution 



for r > _R is 



r\T) = h{nr). 



In the Coulomb limit (22 1 reduces to f°^^(r) ex r 



(22) 

Now let us turn our attention to the plasma region, 
r < R. Here the situation is more complicated since the 
dielectric function depends on the radial distance from 
the trap center. Writing the radial function as 



f{x,uj) = X g{x,uj), 
leads to the following equation for g{x,iLi), 
d 



(23) 



x-rT- [e{x,uj)g'{x,uj) 



2{l+l)e{x,uj)g'{x,uo) (24) 
-[x — £e'{x, ll>)] g{x, w) — 0. 



We now perform another change of variables from x to 
z via x^/xl = z with xl = -2172ei = 2(^1(0) - f^^), 
accompanied by g{x,u!) — ?► g{z,rt). Using the explicit 
result (17) for the dielectric function, Eq. (24) turns into 



a hypergeometric differential equation for g{z, 51) 

z(l - z)g"{z, n) + [i + 3/2 -{£ + 5/2)z] g'{z, n) (25) 

— ^.9(z,r!) = 0, 



which has the general solution 

• — Si ai + Si 



~g{z,n)^A2Fi 



ai;z ] + 



Bz'^-°" 2F1 



Pi -Si (3i + , 



-;(3i;z 



(26) 



around z = 0. Here A and B are arbitrary constants and 
the parameters of the hypergeometric function 2F1 are 



ai = e+-, I3i = --i, Si 



£{£ + !) + -+2Q^. 



From Eqs. ( 23|26 ) we obtain two independent solutions 
of Eq. ( 20 ) inside the plasma. 



/W(x,L.)=xSi^l 



f ai- 5i ai + 5i 



-]oii; 



f^^\x,^)=x-('+^\F, 



Pe-Si Pi + Si 



;A;- 



Only /'^•'(a;) is finite at the origin and thus constitutes 
the correct solution for r < R. 

Collecting the previous results, the solution inside the 
plasma is given by 



An, N I rp fai-Se ai + Se kV^ 



(27) 



The hypergeometric function describes how the per- 
turbed potential is modified for a Yukawa plasma when 
compared with the solution for Coulomb interaction, 
where /'"(r) = r^. 

Since the normal modes explicitly depend on the 
eigenfrequencies their discussion will be postponed to 
SeciniDl 



C. Eigenfrequencies 

1. Existence of an upper bound for the eigenfrequencies 

Before we explicitly discuss the eigenfrequencies we in- 
spect Eq. ^ more closely. Following Ref. [57] we mul- 
tiply by 0* and integrate over volume. Using Gauss's 
theorem we obtain 



n^\(j)i\^ + eiKr,uj)\V(j)i\-' ) dr = 



(28) 



The first integral is always positive which implies that 
there must be a region w here e(Kr, w) < 0, i.e. < 



trivial solutions. Since we used an expansion in spher- 
ical harmonics for 01, Eq. (19a I reduces to the conti- 



considered here the maximum plasma frequency is at the 
center, max[r2p(x)] — fip(O). 

Explicit results for the undetermined eigenfrequencies 
are now found by requiring Eqs. ( 19a 19b I to yield non- 



ri^ < max[r2^(x)], see Eq. (17). For the density profile nuity of /(r,a;) and Eq. (|19b|) requires the continuity of 



e(Kr, w)/'(r, w) across the plasma boundary for any given 
i. The necessary condition for a non-trivial solution is the 
vanishing of the determinant. 



^m-c 






2F1 






I 



ae + Se ^i ^'^ n 

1,^^— + l;a, + l;— j =0, 

(29) 



where we used the property 



d2Fi{a,b;c;z) 
dz 



ah 



i^i(a + l,fe+l;c+l;z). 



Note that Eq. (29) only involves Vl^ and ^. The eigen- 
frequency uj can easily be extracted from il^ by solving 
il^ — uj{uj + iv) for w, which yields the same relation 
as for the normal modes in the discrete A'^-particle sys- 
tem (damped harmonic oscillator) [28]. In the absence 
of damping we have O^ = cj^/wq, i.e. VI reduces to the 
eigenfrequency scaled by the trap frequency. Before we 
proceed to the solutions of Eq. ( 29 ) let us discuss some 
properties of the hypergeometric series |26j 



iFi{a,b;c;z) 



y^ (a)fc(b)fc z'' 

fe=0 



= 1 



(c)fe k\ 

ab z a{a + 1)6(6 + 1) z^ 
VTI ^ c(c+l) 2\ 



(30) 



where {a)k = r(a -I- fc)/r(a) denotes the Pochhammer 
symbol. Its convergence is assured for \z\ < 1 if c is not 
a negative integer and for \z\ = 1 if ^{c — a — b) > 0. In 
our case the condition 3(?(c — a — b) > is not satisfied as 
c — a — 6 = (or —1 for the derivative). 

We now apply this to Eqs. ( 27|29 ). The convergence 
condition \z\ = \K^r'^ /x1\ < \^/x^< 1 is closely con- 
nected to the plasma frequency at the boundary. 



KiO = s 



e 



e>o. 



(31) 



since it is fulfilled for all il^ < nl{£,). 

To better understand the nature of this maximum fre- 
quency we recall the Coulomb limit. In this case the 
density is r independent and there exis ts a unique plasma 
frequency, ftp = -s/S, defined by Eqs. (8|9), at which the 
dielectric function (17) vanishes. This is in contrast to 



Yukawa interaction due to the inhomogeneous density 
profile. However, the frequency (31) is the plasma fre- 
quency at r = i? (a; = ^), i.e. it is precisely the one 
for which the dielectric function vanishes at the plasma 
boundary, e(^, J7p(^)) = 0. 



For any fl^ in the interval f7^(0) > rj^ > rj^(^) there 
exists a point Xs inside the plasma at which 51^ = ^^(xs) 
and consequently e{xs,^) — 0. At this point the lo- 
cal plasma frequency is in resonance with the mode fre- 



quency and the differential equation (20) has a singular 
point. In these cases the solutions of (20) can be singu- 



lar and are associated with a continuous spectrum and, 
possibly, damped quasi-modes, see e.g. [37|[1S]. The ap- 
propriate approach for uncovering the quasi-modes is a 
Laplace transform of Eqs. (|4|. Here we restrict ourselves 
to regular normal mode solutions with Q^ < Op(^). 



2. Coulomb limit and frequency degeneracy lifting of bulk 
modes in a Yukawa plasma 

Having obtained these general properties we now ex- 
plicitly determine the eigenfrequencies. Let us first con- 
sider the limit ^ ^ 1. Performing a series expansion of 
Eq. (29) ioT £ ^ (details can be found in Appendix [b|) 
we find 

2^+1 ^ V 8£3 + 12^2 _ 2^ - 3 , 



nm 



e 



(32) 



For Coulomb interaction {S, = 0) we get nj = 3^/(2^+1). 
This is the well known result for surface oscillations of a 
homogeneous plasma sphere |30j. Furthermore it is eas- 
ily verified that J7i = 1 is a solution of Eq. (29) for any 
^. This mode describes the center of mass oscillation 
(dipole or sloshing mode) and is independent of the par- 
ticle number and the screening parameter |31j . 

Another solution in the Coulomb limit is given by the 
Coulomb plasma frequency flp. This can easily be seen 
from Eq. ([6]) with k = 0. Since the dielectric function 
inside the plasma vanishes for fl = Up, Eq. (|6| is satisfied 
for any 0']" that satisfies the boundary conditions (19). 



This implies a high degeneracy since these requirements 
can be met by an infinite number of modes with arbitrary 
£, m. 

For Yukawa interaction with a finite ^ this degeneracy 
is lifted and a series expansion yields (see Appendix [B|) 

f},2(^)~3 + c,e' + ..., C«l, (33) 



where cq « 0.85031, ci « 0.98624 and C2 w 0.99992. 
From this series of expansion coefScients we see that the 
frequency approaches flp(^) ~ 3 + ^^ + . . . as ^ increases. 
For fi^ = ^p(0 the hypergeometric series does not 
converge in general at \z\ = 1. However, a closer in- 
spection of the series representation ([30]) reveals that a 
well behaved solution can be obtained if a or 6 is a neg- 
ative integer — (n — 1). In these cases the series is sim- 
ply a polynomial of order n — 1, where n G N^ . Since 
b — {ai+Si)/2 > we require a — (a^ — (5f)/2 = — (n — 1). 
Solving this equation for ft^ yields 



nl,{CnT)^i2n + l)in-l) + i2n-l)£. 



(34) 



Keeping in mind that at this point ^p{£,) = 3 + ^^/(1 +^) 
we find 






1 



Cnl + VCni iCne + 4) 



C„« = (2n-l)(n + £)-4, 



(35) 



where i > 3 for n — I. Even though this choice guar- 
antees a well behaved solution, it can be shown that 
Eqs. ( 34|35 l do not solve the eigenvalue equation (29 1. 
Nevertheless, a numerical evaluation shows that solutions 
do exist for ^ > ^^"* which closely approach ^p{£,) as 
? — >■ Ci"*- This issue will be dealt with in more detail in 
SeclmDl 



3. Macroscopic/strong screening limit, ^ ^ 1 

Let us now discuss the limit ^ — > oo. Using fl^ < flp{^) 
one can show that I'^^/a;^] — > 1. As before, the con- 
vergence problem at this point can be circumvented by 
choosing the parameters of the hypergeometric function 
such that its series terminates at a finite order. Requiring 
a = [oLi — 5i)/2 = —n (n S N) and solving the equation 
for H,^ yields the solution for the eigenfrequencies in the 
limit ^ — > oo. 



^— >-oo 



n 



2n'^ + i2e + 3)n + £. (36) 



The reason we chose the same index n as in the previous 
case will become clear shortly. 



It is shown in Appendix p5l that ( 36 ) actually solves 



Eq. (|29|). Further, the lowest order correction for finite ^ 



is found as 



nt 



e 



+ .. 



?»i, 



(37) 



where the coefficients are given by 



dni = (2n + 1 + 3/2) [(4n^ + 12^^ + 3n - 9)n 

-t- 2^(471^ + 8^2 + 2tn{ji +!)+(. -I)]. (38) 




FIG. 2: Eigenfrequencies and their dependence on ^ for var- 
ious modes {n,£). Also shown (by the arrows) are the limits 
for ^ = and ^ — >■ oo. The crosses denote the parameters at 
which new modes appear to the right of ^Jj"' . 



4- Eigenfrequencies for arbitrary ^ 



For arbitrary values of ^ we solved Eq. ( 29 ) numer- 
ically. The results are shown in Fig. [2J In the weak 
screening limit, ^ <C 1, the known Coulomb limit is recov- 
ered, where fi^ = 3£/{2£ + 1) (surface modes) or fl^ — 3 
(bulk modes) [30] ■ Except for the center of mass mode 
all mode frequencies increase with ^ and saturate in the 
limit ^ — > oo. We find numerically that the eigenmodes 
with mode number n > 1 (£ > 3 for n = 1) at ^ = cxd 
approach ^p{^) as ^ is decreased and cease to exist for 
^ < Cn7*- Thus, the chosen indices are the same and the 



modes evolve continuously from C" to £_ 



oo. 



Let us briefly summarize the findings of this section. 
The main result is that screening lifts the degeneracy of 
the bulk modes of the Coulomb system and the appear- 
ance of a new mode number n. For a given value of f 
the number of allowed modes is restricted. Modes with 
n = {£ > 1) and n = 1 with £ = 0, 1,2 exist for all 



^ > whereas our numerical solution of Eq. ( 29 1 indi- 
cates that all other modes exist only for ^ > ^^- The 
(n, £) = (1, 3) mode is a special case with C"* = 0. 

Having found the eigenfrequencies we can now come 
back to the discussion of the shape of the normal modes. 



D. Explicit results for the normal modes 

1. Coulomb limit 

The eigenmodes in the Coulomb limit are well known, 
see e.g. J^. The surface modes with Clj = 3^/(2^+1) are 
given by ^Tl^ r-('^+i)r™(6', </?) and 0i" - r%"'{e,ip), 
ci. Eqs. ( 22|27 ), while the bulk modes oscillate at the 
plasma frequency f2p = v3- The potential eigenfunc- 



tions in the latter case are undefined inside the plasma, 
which can easily be seen from Eq. m. Since the dielec- 
tric function in the Coulomb limit is constant for r < R, 
Eq. (k)) is satisfied for any (f>"'^ if J7 = Qp. The potential 
pertu rbat ion is only restricted by the boundary condi- 
tions ([l9|. It follows from Eqs. ^ and A0°"* = that 
07' = 22|. This further imphes that /'"(i?) 
Eq. (pal).' 



0, see 



In the following wc will discuss the eigenmodes for 
Yukawa interaction and point out the similarities and 
differences when compared with the Coulomb limit. 



p: 0.5 
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FIG. 3: Radial eigenfunctions in the limit ^ — > oo for various 
modes as indicated in the figure. 



2. General remarks for ^ > 



For ^ > the radial eigenfunctions (27) inside the 
plasma explicitly depend on fl^ , which was shown to have 
several solutions for a given £. Thus, in addition to the 
angular mode numbers m and £, there is a radial mode 
number n which determines the structure of the radial 
eigenf unction. The radial function is given by Eq. (27) 



and the corresponding eigenfrequency is determined by 
Eq. (l29l). 



3. Eigenmodes for Yukawa interaction 

Let us begin with the limit ^ — ^ oo. It was already 
shown in Sec. Ill C| that in this case the hypergeometric 
function reduces to a polynomial of order n in z — r^ /R^. 
From Eqs. (27|36) we thus get the radial eigenfunctions 



rir,t.^,)^r' 



k\ {a,)k \r) 



(39) 



fc=0 

OC /p(^+l/2.0) 



2r2 



where the Pn ' {x) are Jacobi polynomials. The 
n = surface modes are particularly simple. Here, the 



sum in Eq. (39) is just a constant and the eigenmodes 



inside the plasma are the same as in the Coulomb limit. 
There are no radial nodes. For finite ^ we find that the 
solutions differ only slightly from the results at i^ = or 
(XI. These modes are the natural generalization of the 
surface modes to Yukawa interaction and there are no 
qualitative changes compared to a Coulomb system. In 
particular, the dipole modes with {n,£) = (0, 1) describe 
the three center of mass oscillations and have eigenfunc- 
tions (/)'" ^ rY{"{9, ip), independent of ^. 

The origin of the modes with n > can be traced back 
to the bulk modes of the Coulomb system, see Fig. [2] 
While for 1^ = the eigenmodes are not entirely specified, 
for ^ > their form is determined by Eq. (pTJ), together 
with ([29]) for the eigenfrequency. In the limit ^ — j' 00 
the potential perturbations are given by Eq. (39), see 
Fig. [3J n is the number of radial nodes. With increasing 



£ the nodes and extrema are shifted towards the cluster 
boundary. 

The behavior for finite ^ is shown in Fig. |4j The n = 1 
modes with £ = 0, 1, 2 extend up to ^ = 0, while all other 
modes with n > 1 exist only for ^ > ^™*. 

Let us discuss the former case first. As we move from 
the f = 00 limit towards <^ = 0, the single node ap- 
proaches r = R, see Fig. |4^,b. The radial eigenfunctions 
at ^ = 0+ read 



r(r,0 = rSFi 



V 2 ' 2 '"^'3-2c£ 



where £ = 0,1,2 and Sg must be evaluated at il^ — 3. 
Compare also with the series expansion performed in Ap- 
pendix pi This form is in accordance with /'"(i?) = 
for ^ = (the coefficients cg were chosen to ensure this). 
Thus, these modes also exist in a Coulomb system, where 
they are among the bulk modes with Q = flp. The main 
difference between Coulomb and Yukawa interaction at 
this point is that in the former case the mode form is not 
specified by Eq. ^. For C > this equation constitutes 
an additional restriction on the mode form and selects 
the eligible modes for ^ = 0+ from the large number of 
modes at ^ = 0. 

In the case of all other modes with n > 1 the outermost 
node comes arbitrarily close to r = i? as ^ decreases, but 
disappears at CS'*' ^^^ ^^&- Ppjd- This goes along with 
a very strong increase of dP'^{r)/dr at the boundary, 
which makes the numerical solution of Eq. ( 29 ) increas- 
ingly difficult. If 17 = rjp(^™*) were a proper solution, 
the boundary conditions \19\ and Eq. (21 ) would require 
/°"' == 0, and hence /'"(i?) = 0, just like for ^ = and 
il = rip. However, the potential inside the plasma is not 
undefined but determined by the solutions of Eq. (|6|. 
Analogously to the case ^ = 00, Eq. (27) reduces to a 
polynomial at ^p{£,nt) ^^^ cannot satisfy the boundary 
conditions since f™{R) ^ 0. Our numerical solutions of 
Eq. ( 29 ) indicate that these modes exist only for ^ > ^J;^'* ■ 
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FIG. 4: Radial eigenfunctions fnt{r) for various ^ and mode 
numbers {n,£). 
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FIG. 5: Radial fluid velocity for the breathing mode. 



Jf.. Breathing mode 

The lowest monopole mode has the index (n, () — (1, 0) 
and is worth a more detailed discussion. For ^ — > oo the 
fluid velocity satisfies v™ ~ re^ since 0']" is quadratic in 
r, see Eq. (5b). This corresponds to a uniform breath- 
ing oscillation of the plasma. In the Coulomb limit this 



mode is among the bulk modes. It was shown in Ref. [35] 
that a system of N harmonically confined particles with 
Yukawa interaction does not support a universal uniform 
breathing mode. We find a similar behavior in our fluid 
approach since the uniformity condition is only fulfilled 
for ^ — > cx). For finite ^ the uniformity condition is not 
satisfied, see Fig. [5j Further, it is surprising that it is 
not recovered as ^ — >■ 0. However, this is not contradic- 
tory since this solution is just one among many in the 
Coulomb limit. Here, the breathing mode is given by 
/'"(r) ^ (r^ — i?^), which also satisfies the boundary 
condition /'"(i?) =0 Hg. 

The aforementioned breathing mode has been studied 
by Sheridan |23j under the assumption of a homogeneous 



density and a uniform displacement. We find that his an- 
alytical result for ^ (^) agrees with our numerical results 
to within 0.2% for any ^. Note, however, that the equa- 
tion ^ is determined from in [53] is different from Eq. ( 10 1, 



since ^ used in [23 involves the radius for a homogeneous 
sphere. This mode was found to be the dominant ex- 
citation during spherical crystal formation after a rapid 
temperature quench ;16j . 



IV. CONCLUSION 



Summarizing our results, we have investigated the 
ground state and the normal modes of a harmonically 
confined Yukawa plasma in a fluid approach. The results 
of TT] were extended to a time-dependent theory and 
simple formulas for the plasma radius, total energy and 
density moments were obtained. The density moments 
could be used to infer the dimensionless plasma parame- 
ter ^ = kR, the ratio of the plasma radius to the screening 
length, from experimental data. For typical dusty plasma 
experiments with na « 0.5 ... 1 and N ^ 50 . . . 1000 we 
expect ^ « 1.5 ... 6 [15i. 

Further, the fluid equations were linearized and solved 
for the normal modes. Compared with previous results 
for Coulomb interaction, we found a new class of modes 
with a radial mode number n that determines the num- 
ber of radial nodes in the potential eigenfunctions. The 
eigenfrequencies were found to depend only on ^. The 
degeneracy of the bulk modes in the Coulomb limit was 
shown to be lifted for Yukawa interaction and series ex- 
pansions for the eigenfrequencies were derived for limit- 
ing cases. For experimentally relevant parameters they 
must be determined numerically, though. 

The fluid theory should be applicable for large clusters 
with N > 100. In [IT] it was shown that the agreement 
between the continuum theory and the exact A'^-particle 
ground state was good at low screening, but deviations 
increased for larger Ka. Applying the local density ap- 
proximation, the authors were able reduce the deviations 
for strong screening [T5]. The same behavior was ob- 
served in [16j , where the frequency of the monopole mode 
excited in a Langevin dynamics simulation was compared 
with the theory of |23j . Thus, similar behavior is ex- 
pected for the results presented here. Additionally, one 
has to bear in mind that the fluid equations correspond to 
a mean-field description and neglect correlation effects. 
They are thus not able to describe crystallization and 
shell structure formation. The same applies to the local 
density approximation used in [T5] . This question was re- 
cently analyzed in [33j [34] . Only a comparison with first 
principle simulation data will show the true applicability 
limits of the fluid approach. This analysis is subject of 
ongoing work. 
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Expanding the squared eigenfrequency as Q,^ (0 ~ "^ ' 
bS^ + c£,^ we find 



Appendix A: Calculation of (jjoir) 



Carrying out the angle integration in Eq. (3b) wc ob- 
tain [ri] 



00 (r) = 27r— / dr'no{r')r' 
Kr Jo 



g-K|r-r'| _ g-K(r+r') 



The density involves a constant term and a term ^ r 
Thus we consider the following integrals for n — I, 3, 






dr' {r'Y 



g-K|r-r'| _ g-K(r+r') 



li r > R we have |r — r'| = r — r' for the entire integration 
and 



■ / dr' (/)" 
Jo 



n = 1, 



^ (Kcosh(o-sinh(e)]«^ 

\C(^' + 6) cosh(^) - 3{e + 2) sinh(0, n = 3. 

For r < R the integrals must be solved independently in 
regions where r' < r and r' > r, 



I<^I>{i^nT)+ / d/(/) 



Z> (C ^ Kr) + 2 sinh(Kr) / dr' (r')"e 



g-K(r'-r) _ g-K(r+r') 



^' { ,^'\^^^^'^ 



e(0«i- 



36 \ 



e 



3ac - 36^ - a2 



e 



Further, we obtain x^ sa 6 — 2a and —^k'^{^)/ki{S,) ~ 
(^ + 1) + CV(2^ - 1) for £ ^ or -^fc^(6Ao(0 = 1 + 
^. The series expansion for the hypergeometric function, 
Eq. (30), then yields, for a ^ 3, 



2^1 



ae -Si ae + Se ^^ 



a — 3 / 6 



e 



It is sufficient to approximate the other hypergeomet- 
ric function in Eq. (29) by 1. Comparing terms of or- 



der £,'^ , S,^ , i,"^ , we find the coefficients a,b,c as given in 
Eq. ([32]). 

If a = 3 we must keep in mind that fl'^{S,) < ^piS,)- The 
ansatz n'^{^) « 3 + 6.^ yields ^V^^s ~ -'?/(2^)- However, 
this choice turns out to be inadequate since it leaves a 
constant term in Eq. (29). Next, we use ri^(^) « 3H-c^^, 
which results in C'^/x^w 1/(3 — 2c). Then, the lowest 
order term in Eq. ( 29 1 vanishes if 



ai- Si ai + S/, 

2-^1 1 — ^ — , — 7, — ;a^; 



1 

3 -2c 



^2=3 



This equation can be solved numerically for c if £ = 0, 1, 2 
and yields the coefficients given below Eq. ( 33 ) . For ^ > 3 



we find no solution. The coefficients are in accordance 
with the condition 57^ (^) < f^p(C) ~ 3 + ^^ as they satisfy 
c<l. 



The remaining integral can easily be solved. Collecting 
the results and using the explicit result for no(r) [Eq. ^] 
we obtain Eq. (14). Alternatively, the potential inside 



the plasma may be directly calculated from Eqs. (7,10) 

ofUB. 



Appendix B: Series expansion for f^^(C) 

The difficulty in finding an expansion for i^^(^) arises 
from the hypergeometric function 2^^1(0, b; c; z) since a, b 
and z are all functions of ^. In the following we will 
separately discuss two different limits. 



2. Macroscopic/strong screening limit, ^ 3> 1 

In this limit we seek an expansion in terms of y 
^~^ <^ 1. The ansatz J^^(y) ^ d + by + cy"^ yields 



e{y-') 



1 -a(d-2)+b y-2 

-— H r^ , ^wl-2y, 

ay a^ a;^ 



and —k'^{y ^)/ki{y ^)/y fa y ^ + 1. The expansion for 
the hypergeometric function is found from |26j 
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iFi (a, 6; a + 5 — m] z) 



(to — 1)! r(a + 6 — in) 



m — 1 



T{a)T{b) 



(1-z)-™^ 



(a-m)fe(6-TO)fe(l-z)'= (-l)"r(a + 5-TO) 



fc=0 



fc!(l - m)*; 



r(a — TO,)r(5 — TTi) 



E 

fe=0 



(Q)fc(&)fc 
fc!(A; + TO)! 



(- ln(l - z)+ i;{k + 1) + V'lfc + TO + 1) - V'la + k) - tP{b + fc))(l - z)"" , 



valid for jz — 1| < 1 and to e N+. ijj{x) denotes the 
Digamma function. The expansion for m = Ois the same, 
but without the first (finite) sum. 

The following calculation is based on a = flf^^ ^ , see 
Eq, ^ 



36| , and will show that this choice solves Eq. ( 29 1 . 
For the first hypergeometric function in Eq. ( 29 1 we have 
m — 0, a ^ —n — q{y) and 6 « n + ae + q\y), where 
l{y) "^ 1- The main contribution in the sum arises from 
■0(0 + fc) « '4'{—n + k — q) « q~^ (for k < n), which 
yields a constant term when combined with the prefactor 
r-i(a) !=ir-\^n~q) 
term is then found as 



(— l)"+^n!(7. The lowest order 



2Fi{- n - q,n + ai + q;ai;l - 2y) 



(-1)" 



n\. 



Similarly, we obtain for the hypergeometric function with 
shifted parameters and n = 1,2,..., 

2Fi{ - n + I - q,n + I + at + q] at + l;l - 2y) Ki 



ai 



(-1)" , 
{ai)n 



2xn{n + ai) 



X = An + 2ai 



where we used the above expansion for m = 1 and q{y) « 
by/X- In the following the notation 2-^1 [+1] (2-Fi[+0]) 



will be used to denote the (un) shifted hypergeome tric 
function. Comparing terms 0{y~^) we find that Eq. (29) 
is satisfied if 6 = 0. For n = the previous equation is 
not valid and 2-Fi[+l] kj (2y)^^. In this case the leading 
order term in Eq. (29 1 is 0{y^'^). It vanishes due to its 
prefactor ^ {£ — a), since a = ^ for n = 0. The 0{y~^) 
term vanishes for 6 = 0. 

In order to calculate the lowest order correction to 
r2^(j/) for n — 1,2,... we need to evaluate the hyper- 
geometric function up to first order in y. We find, using 
q{y) « cyVx, 



2F1 [+0] « ^^n! [1 - 2n{n + ai)y] , 2^1 [+1] ~ a^ 



(-1)" , 
{ai)n 



{0lt)n 

-1+ ( (n- l)(n + af + 1) + 



2xn{n + ai) 



The coefficient c can now be determined by choosing it 
such that terms 0{y^) vanish, yielding Eqs. ( 37|38 ) for 
n> Q. If n = the same procedure leads to c = —2t{(. — 
l)(^ + 3/2), i.e. Eq. (|38| also holds for n = (here £ 7^ 0). 
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